/*
 * fdmconstants.h
 *
 *  Created on: Apr 11, 2010
 *      Author: nish
 */

#ifndef FDM1DUTILS_H_
#define FDM1DUTILS_H_

#include "fdmconstants.h"

double length = 5e-10;
//double a = 5e-12;
double a = 0;
double t0 = 0;
int npoints = 10;

double hartree_Const = 0;

void printConstants();
void initConstants();
void loadConstantsFromFile();
gsl_matrix* formHMatrix();
void printHMatrix(gsl_matrix* H);
void saveMatrix(gsl_matrix* H, char* name);
void saveVector(gsl_vector* V, char* name);
void solveEig(gsl_matrix* H,gsl_vector* eig_val, gsl_matrix* eig_vect);
void printEig(gsl_vector* eig_val, gsl_matrix* eig_vect);
void makePlotFiles(gsl_vector* eig_val, gsl_matrix* eig_vect, char* name_vect, char* name_val, char* name_prob);
gsl_vector* getProbabDensityFromEigVect(gsl_matrix* eig_vect, int index);
void printVector(gsl_vector* vect, char* name);

/* Hydrogen 1D atom spherical */
gsl_matrix* formHydrogenHMatrix(gsl_vector* pot);
gsl_vector* formPotHydrogen1D();

/********** HE ATOM *************/
gsl_matrix* formHeliumHMatrix(gsl_vector* pot, gsl_vector* probabDensity);
gsl_vector* formPotHelium1D();
double getHartreePotential(gsl_vector* probabDensity, int index);
gsl_vector* formHarteePotential(gsl_vector* probabDensity);
gsl_vector* formInitialHartreePotential();
double getErrorFromHartreePotentials(gsl_vector* oldHartree, gsl_vector* newHartree);

/********************************/



/*** Hydrogen 1D methods ***/



/**
 * 1. Form the Hmatrix -> involves the 1D spherical coordinates ->
 * include the nuclear potential and spehrical harmonics.
 * Ef(r) = (-hbar^2 * (d^2/dr^2) * 1/2m - q^2/4*pi*eps*r + (l(l+1))*hbar^2/(2mr^2)) f(r)
 * for the nuclear potential and the sperical harmonics term, take r as starting from
 * a to 100*a, where a = 0.1A and 0.05A
 *
 */

gsl_matrix* formHydrogenHMatrix(gsl_vector* pot){
	gsl_matrix* H = gsl_matrix_alloc(npoints,npoints);
	int i = 0;
	int j = 0;
	double nuclr_pot = 0;
	double sphr_hrmonic = 0;
	double oneDterm = 0;
	double a_incr = a;
	int l = 0; // s-orbital for p, d and above orbitals -> use l =1,l=2 etc


	for(i=0;i<npoints;i++){
		a_incr = a;
		if(i== npoints){
			break;
		}

		for(j=0;j<npoints;j++){
			if(j== npoints){
				break;
			}

			if(i==j){
				a_incr+=a*j;
				oneDterm = 2*t0;
				nuclr_pot = gsl_vector_get(pot,j);
				sphr_hrmonic = l*(l+1)*pow(hbar,2)/(2*e_mass*pow(a_incr,2)*e_charge);
			}else if(abs(i-j) == 1){
				oneDterm = -1*t0;
				nuclr_pot = 0;
				sphr_hrmonic = 0;
			}else{
				oneDterm = 0;
				nuclr_pot = 0;
				sphr_hrmonic = 0;
			}

//			printf("The varying:%d, %.5g, %.5g, %.5g\n",j,nuclr_pot,sphr_hrmonic,a_incr);
//			printf("The varying:%d, \n",j);

			gsl_matrix_set(H,i,j,(oneDterm+nuclr_pot+sphr_hrmonic));
		}
	}
	return H;
}

gsl_vector* formPotHydrogen1D(){
	gsl_vector* pot = gsl_vector_alloc(npoints);
	int i=0;
	double a_incr = a;
	double nuclr_pot;
	FILE* f = fopen("pot_data-H1d.dat","w");

	for(i=0;i<npoints;i++){
//		if(i==0){
//			nuclr_pot = -1*1e-19;
//		}else{
			nuclr_pot = -1*(e_charge)/(4*M_PI*eps*a_incr);
//		}
		gsl_vector_set(pot,i,nuclr_pot);
		a_incr+=a;
		fprintf(f,"%d\t%.5g\n",i,nuclr_pot);
	}

	fclose(f);

	return pot;
}

/******************************/

/************************************** He Atom with Hartree potential *********************************/


gsl_matrix* formHeliumHMatrix(gsl_vector* pot, gsl_vector* hartreePot){

	gsl_matrix* H = gsl_matrix_alloc(npoints,npoints);
	int i = 0;
	int j = 0;
	double nuclr_pot = 0;
	double sphr_hrmonic = 0;
	double oneDterm = 0;
	double hartree = 0;
	double a_incr = a;
	int l = 0; // s-orbital for p, d and above orbitals -> use l =1,l=2 etc


	for(i=0;i<npoints;i++){
		a_incr = a;
		if(i== npoints){
			break;
		}

		for(j=0;j<npoints;j++){
			if(j== npoints){
				break;
			}

			if(i==j){
				a_incr+=a*j;
				oneDterm = 2*t0;
				nuclr_pot = gsl_vector_get(pot,j);
				sphr_hrmonic = l*(l+1)*pow(hbar,2)/(2*e_mass*pow(a_incr,2)*e_charge);

				/*
				 * Call the hartree function here
				 */
//				hartree = getHartreePotential(probabDensity,i);
				hartree = gsl_vector_get(hartreePot,i);

			}else if(abs(i-j) == 1){
				oneDterm = -1*t0;
				nuclr_pot = 0;
				sphr_hrmonic = 0;
				hartree = 0;
			}else{
				oneDterm = 0;
				nuclr_pot = 0;
				sphr_hrmonic = 0;
				hartree = 0;
			}

//			printf("The varying:%d, %.5g, %.5g, %.5g\n",j,nuclr_pot,sphr_hrmonic,a_incr);
//			printf("The varying:%d, \n",j);

			gsl_matrix_set(H,i,j,(oneDterm+nuclr_pot+sphr_hrmonic+hartree));
		}
	}
	return H;
}

double getHartreePotential(gsl_vector* probabDensity, int index){
	double hartree = 0;
	int i=0;

	for(i=0;i<npoints;i++){

		if(i==index){
			continue;
		}
		hartree += (gsl_vector_get(probabDensity,i)/(abs(index-i) * a));
	}


	return hartree;
}

/*
 * Forms the Hartree potential vector for a given probability density
 */
gsl_vector* formHarteePotential(gsl_vector* probabDensity){
	gsl_vector* Hartree = gsl_vector_alloc(npoints);
	double h = 0;
	int i,j;


	for(i=0;i<npoints;i++){
		h=0;

		for(j=0;j<npoints;j++){
			if(j==i){
				continue;
			}
			h+=(gsl_vector_get(probabDensity,j)/(abs(i-j) * a));
		}
		gsl_vector_set(Hartree,i,h*hartree_Const);
	}
	return Hartree;
}

/*
 * Forms a dummy Hartree potential for the first iteration
 */
gsl_vector* formInitialHartreePotential(){
	gsl_vector* Hartree = gsl_vector_alloc(npoints);
	int i=0;
	for(i=0;i<npoints;i++){
		gsl_vector_set(Hartree,i,0);
	}

	return Hartree;
}


/*
 * gets the sum abs difference between all the elements of the old and new hartree potentials
 * and averages them by dividing with npoints. This is the final error/difference
 * in the old and new hartree functions.
 */
double getErrorFromHartreePotentials(gsl_vector* oldHartree, gsl_vector* newHartree){
	double diff = 0;
	double ret = 0;
	int i=0;

	for(i=0;i<npoints;i++){
		printf(">>%.5g\t%.5g\n" ,gsl_vector_get(oldHartree,i),gsl_vector_get(newHartree,i));
		diff = diff + abs(gsl_vector_get(oldHartree,i) - gsl_vector_get(newHartree,i));
	}
	ret = diff/npoints;

	return diff;
}



gsl_vector* formPotHelium1D(){
	gsl_vector* pot = gsl_vector_alloc(npoints);
	int i=0;
	double a_incr = a;
	double nuclr_pot;
	FILE* f = fopen("pot_data-He1d.dat","w");

	for(i=0;i<npoints;i++){
		nuclr_pot = -1*(2*e_charge)/(4*M_PI*eps*a_incr);   //we are not doing e_charge^2 as we are considering all the energies in eV.
		gsl_vector_set(pot,i,nuclr_pot);
		a_incr+=a;
		fprintf(f,"%d\t%.5g\n",i,nuclr_pot);
	}

	fclose(f);
	return pot;
}


/****************************************************************************************************/
void printConstants(){
	printf("here are the constants:t0=%.5g, npoints=%d, a=%8.8g\n",t0,npoints,a);
	return;
}


void initConstants(){
	loadConstantsFromFile();

	a = length/npoints;
	t0 = pow(hbar,2)/(2*e_mass*pow(a,2)*e_charge);

	/*
	 * Hartree const = e*e/(4*pi*eps0)
	 */
	hartree_Const = e_charge/(4*pi*eps);
	return;
}



void loadConstantsFromFile(){
	//read a properties file - for npoints and length config parameters.
	FILE* f = fopen("config.prop","r");
	char line[128];
	char * tokens;

	if(f != NULL){
		while(fgets(line,sizeof(line),f) != NULL){
//			printf("%s",line);
			tokens = (char*)strtok(line,"=");
//			printf(">>%s\n",tokens);

			if(strcmp("length",tokens)==0){
				tokens=(char*)strtok(NULL,"=");
//				printf("length is %.5g",atof(tokens));
				length=atof(tokens);
			}else if(strcmp("npoints",tokens)==0){
				tokens=(char*)strtok(NULL,"=");
//				printf("\nnpoints is %d",atoi(tokens));
				npoints = atoi(tokens);
			}
		}
	}

	fclose(f);

}


void printVector(gsl_vector* vect, char* name){
	printf("----Vector - %s---------\n",name);
	int i=0;
	for(i=0;i<npoints;i++){
		printf("%f\n",gsl_vector_get(vect,i));
	}
	printf("-------------\n");
}

gsl_matrix* formHMatrix(){
	gsl_matrix* htemp = gsl_matrix_alloc(npoints,npoints);
	int i,j;
	gsl_matrix_set_zero(htemp);

	for(i=0;i<npoints;i++){
		for(j=0;j<npoints;j++){
			if(i==j){
				gsl_matrix_set(htemp,i,j,2*t0);
			}else if(abs(i-j) == 1){
				gsl_matrix_set(htemp,i,j,-1*t0);
			}
		}
	}
	return htemp;
}

void printHMatrix(gsl_matrix* H){
	int i,j;

	printf("in the printHMatrix, %d\n",npoints);
	for(i=0;i<npoints;i++){
		for(j=0;j<npoints;j++){
			printf("%.5g\t",gsl_matrix_get(H,i,j));
		}
		printf("\n");
	}
}

void saveVector(gsl_vector* V, char* file){
	FILE* f = fopen(file,"w");
	int i=0;
	double x = 0;

	for(i=0;i<npoints;i++){
		fprintf(f,"%.5g\t",x);
		fprintf(f,"%f\n",gsl_vector_get(V,i));
		x+=a;
	}
	fclose(f);
}


void saveMatrix(gsl_matrix* H, char* file){
	FILE* f = fopen(file,"w");
	int i,j;

	//gsl_matrix_fprintf(f,H,"%.5g");
	for(i=0;i<npoints;i++){
		for(j=0;j<npoints;j++){
			fprintf(f,"%8.8g\t",gsl_matrix_get(H,i,j));
		}
		fprintf(f,"\n");
	}

	fclose(f);
}


void solveEig(gsl_matrix* H, gsl_vector* eig_val,gsl_matrix* eig_vect){

	//initialize the eigen vect workspace and get solve for eigen values and vectors:
	gsl_eigen_symmv_workspace* workspace = gsl_eigen_symmv_alloc(npoints);
	gsl_eigen_symmv(H,eig_val,eig_vect,workspace);
	gsl_eigen_symmv_free(workspace);
	gsl_eigen_symmv_sort(eig_val,eig_vect,GSL_EIGEN_SORT_VAL_ASC);
//	FILE* f = fopen("eigValH1d.dat","w");
//	gsl_vector_fprintf(f,eig_val,"%.5g");
//	fclose(f);
//	f = fopen("eigVetH1d.dat","w");
//	gsl_matrix_fprintf(f,eig_vect,"%.5g");
//	fclose(f);

	return;
}

void printEig(gsl_vector* eig_val, gsl_matrix* eig_vect){
	int i=0;
	double val;
	gsl_vector_view vect;

	if(eig_val == NULL){
		printf("The eigVal is null");
	}

	for(i=0;i<npoints;i++){
		val = gsl_vector_get(eig_val,i);
		vect = gsl_matrix_column(eig_vect,i);

		printf("Eigen value number:%d -> %.5g\n",i,val);
		gsl_vector_fprintf(stdout,&vect.vector,"%.5g");
	}
		

//	printf("Here is the eigVal size>> %f\n",eig_val->size);

}

/*
 * index is the number of the eigen function for which we want the probability density
 */
gsl_vector* getProbabDensityFromEigVect(gsl_matrix* eig_vect, int index){
	int i=0;
	int j=0;
	gsl_vector* probabDensity = gsl_vector_alloc(npoints);
	double d = 0;

	//index needs to be between 0 and npoints.
	for(i=0;i<npoints;i++){
		d = pow(gsl_matrix_get(eig_vect,i,index),2);
		gsl_vector_set(probabDensity,i,d);
	}

	return probabDensity;
}

void makePlotFiles(gsl_vector* eig_val, gsl_matrix* eig_vect, char* name_vect, char* name_val, char* name_prob){
	int i=0;
	int j=0;
	FILE* f = fopen(name_vect,"w");
	FILE* g = fopen(name_val,"w");
	FILE* k = fopen(name_prob,"w");

	double x =0;
	for(i=0;i<npoints;i++){
		// print the x value
		fprintf(g,"%d\t%.5g\n",i,gsl_vector_get(eig_val,i));
		fprintf(f,"%.5g\t",x);
		fprintf(k,"%.5g\t",x);
		x+=a;
		for(j=0;j<npoints;j++){
			fprintf(f,"%.5g\t",gsl_matrix_get(eig_vect,i,j));
			fprintf(k,"%f\t",pow(gsl_matrix_get(eig_vect,i,j),2));
		}
		fprintf(f,"\n");
		fprintf(k,"\n");
	}

	fclose(g);
	fclose(f);
	fclose(k);

}
#endif /* FDM1DUTILS_H_ */
